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2D  electrostatic  potential  solver  for  Hall  thruster 

simulation 

Justin  W.  Koo* 


This  paper  presents  the  formulation  of  a  2D  axisymmetric  electrostatic  potential  solver 
(2DFV)  for  Hall  thruster  simulation.  It  is  based  on  a  finite  volume  discretization  of  a  cur¬ 
rent  conservation  equation  where  the  electron  current  density  is  described  by  a  Generalized 
Ohm’s  law  description.  Comparison  of  2DFV  to  an  existing  ID  axisymmetric  electrostatic 
potential  solver  is  provided  and  investigation  is  performed  into  anomalous  mobility  correc¬ 
tions  and  Hall  current  calculation.  Details  of  an  extension  of  this  formulation  to  include 
j  x  B  plasma  turbulence  terms  are  also  provided. 


Nomenclature 


j  Current  density,  A/m 2 

ji  Ion  current  density,  A/m 2 

je  Electron  current  density,  A/m 2 

p  Electron  mobility,  m2  /V/s 

B  Magnetic  field,  T 

ne  Plasma  density,  1/m3 

nxe  Neutral  density,  1/m3 

E  Electric  field,  V/m 

p  Plasma  pressure,  Pa 

a  Cross-correlation  terms,  A-T/m2 

h  Normal  vector  (pointing  out  of  cell) 

Te  Electron  temperature,  K 

vth  Electron  thermal  velocity,  m/s 

ven  Electron-neutral  elastic  collision  frequency,  1/s 

Q  Hall  parameter 

mxe  Xenon  mass,  2.15E-25  kg 

me  Electron  mass,  9.11E-31  kg 

qe  Electron  charge,  1.602E-19  C 

kb  Boltzmann  constant,  1.3807E-23  J/K 

dc  Acceleration  channel  width,  0.025  m 


I.  Introduction 

Modern  Hall  thruster  device  simulation  codes,  including  those  by  Fife,1  Hagelaar,2  and  Koo3  use  potential 
solver  formulations  which  rely  on  the  concept  of  the  thermalized  potential  (pioneered  by  Morozov4)  to 
discretize  a  multidimensional  physical  domain  into  a  single  quasi-axial  dimension  demarcated  by  magnetic 
field  lines.  This  formulation  has  proven  to  be  very  successful  but  relies  on  the  strict  assumption  of  isothermal 
electrons  along  magnetic  field  lines  and  is  difficult  to  adapt  to  certain  magnetic  field  geometries.  The 
formulation  of  a  2D  finite  volume  potential  solver  based  on  the  same  concept  of  current  conservation  is 
explored  in  this  paper  and  comparisons  are  provided  to  demonstrate  the  additional  physics  which  can  be 
explored  with  a  true  2D  formulation. 

*  Engineer,  Advatech  Pacific,  Inc.,  Palmdale,  CA,  93550,  and  AIAA  Member. 
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II.  Current  conservation 

Since  2DFV  is  designed  for  use  with  a  hybrid-PIC  code  where  the  timestep  is  much  larger  than  the  time 
necessary  for  the  electric  field  to  equilibrate  to  the  plasma  conditions,  the  simplified  steady- state  version  of 
current  conservation,  shown  in  Eqn.  1,  is  applicable. 

V-J=V-(i;+/e)  =  0  (1) 

When  incorporated  into  hybrid  PIC  codes,  the  ion  current  density  can  be  tracked  directly  via  the 
macroparticles  motion  while  the  electron  current  density  is  derived  from  a  generalized  Ohm’s  Law  for¬ 
mulation. 


III.  Generalized  Ohm’s  law 

The  electron  current  density  is  based  on  the  following  generalized  Ohm’s  Law  formulation: 

je  =  Mie  x5)  +  fineE  +  /iVp  (2) 

Expanding  Eqn.  2,  a  set  of  three  equations  can  be  written  for  the  electron  current  in  the  axial,  radial, 
and  azimuthal  directions  (dropping  the  subscript  e  for  convenience). 

jz  =  p(jeBr  ~  jrB$)  +  nneEz  +  pVzp 
jr  =  l^(jzB$  j$  B  T  /meEr  T 

je  =  p(jrBz  -  jzBr)  +  pneEg  +  pV gp  (3) 

Due  to  the  axisymmetric  geometry  of  the  Hall  thruster,  the  electric  field  and  gradient  of  the  pressure  in 
the  azimuthal  direction  are  identically  zero.  With  this  simplification,  the  azimuthal  electron  current  density 
is  no  longer  a  function  of  any  azimuthal  quantities  and  is  written  as: 


je  =  p{jrBz  -  jzBr) 


(4) 


It  is  now  possible  to  substitute  je  into  the  equations  for  the  axial  and  radial  electron  current  density. 
This  leads  to  the  following  set  of  coupled  linear  equations: 


jz  =  p(p(jrBz  -  jzBr)Br  -  jrBg )  +  pneEz  +  pVzp 

jr  P(jzB$  P(jrBz  j z Br  j B z  J  pneEr  -1-  p'ErP  (5) 

This  equation  system  can  be  solved  to  isolate  the  axial  and  radial  components  of  the  electric  current  density: 


jz  =  Mil  (pneEz  +  pVzp)  +  M12  (pneEr  +  pS7rp) 

jr  =  M21  (pneEz  +  pS7  zp)  +  M22  (pneEr  +  pS7  rp)  (6) 


where, 


Mil  = 

1  +  p2  Bz 

1  +  p2B2 
pBe  +  p2BrBz 

M12  = 

1  +  p2B2 
—fiBe  +  /j?  BrB 

M21  = 

1  +  p2B2 

1  +  p2B2 

M22  = 

1  +  p2B2 
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and, 


Clearly,  in  the  limit  of  Br  =  B  and  Bz 
forms: 


P  = 


Qe 

mev 


(8) 


Bq  =  0,  the  mobility  coefficients  return  to  their  desired  classical 


MuM  MX  =  J  +  ^2^2 

M12M  =  0 
M21M  =  0 

M22M  =  M||  =  M  (9) 


Simplification 

The  following  simplifications  are  used  to  more  easily  manipulate  the  equation  system. 


jz  —  Z3  Ez  +  ZAEr  +  Z-j 
jr  =  B/\  Ez  +  i?4  Er  +  i?5 


Jr  - 

-  jjz  ~r  ~r  j. 

where, 

II 

IN 

enuiz  +  MnMWp  +  M12  MVrp 

^3  = 

/l inline 

II 

Hi2prie 

and, 

enuir  +  fl21^zP  +  P22  p'VrP 

r3  = 

II 

P22  pne 

(10) 


(11) 


(12) 


IV.  Finite- Volume  discretization 


To  formulate  the  numerical  scheme  used  in  2DFV,  integrate  Eqn.  1  over  the  cell  area  as  follows: 


/  V-j<9V  =  /  0dV  =  0 

Jv  Jv 

Using  the  divergence  theorem,  this  equation  becomes: 


(13) 


L 


j  •  hdS  =  0 


(14) 


Eqn.  14  clearly  demonstrates  that  2DFV  simply  balances  the  current  fluxes  into  and  out  of  a  given  cell. 
This  formulation  is  well  suited  to  “noisy”  density  profiles  produced  by  particle  codes  because  it  contains  no 
derivatives  higher  than  first  order  for  the  plasma  pressure  term. 


A.  Computational  stencil 

The  acceleration  channel  of  the  Hall  thruster  (shown  in  Fig.  1)  is  mapped  with  a  cartesian  mesh.  The 
electrostatic  potential  is  stored  at  the  cell-centers  of  this  grid  (nz  x  nr  cells)  while  the  electric  field  values 
(both  axial  and  radial)  are  stored  on  staggered  EW  and  NS  grids  offset  from  the  main  cartesian  mesh  by  Ar/2 
and  Az/2,  respectively.  The  electric  field  values  on  the  EW  grid  (nz+ 1  x  nr  nodes)  are  face-centered  axially 
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and  cell-centered  radially  and  those  on  the  NS  grid  (nz  x  nr  +  1)  are  cell-centered  axially  and  face-centered 
radially.  Discretizing  Eqn.  14  leads  to  the  following  simple  identity: 


jf  A Aeast  ~  jY  AAwest  +  AAnorth  -  j^AAsouth  =  0 


(15) 


Figure  1.  Computational  Stencil  -  (main  figure)  circles  represent  cell-centered  potential  values,  orange  lines 
represents  EW  grid,  blue  lines  represents  NS  grid;  (upper  right  inset)  moving  electric  field  values  from  NS  to 
EW  grid  and  vice  versa;  (lower  right  inset)  purple  lines  represent  cell  normal  directions 

Using  the  simplifications  described  in  Sec.  Ill,  Eqn.  16  can  be  obtained: 


Z%AAeast  +  Z*E§AAeast  +  Z?E%AAeaat  -  Z™ AAwest  -  Z™  E%  AAwest  -  Z™ E™  AAwest 
+R*?  A Anorth  +  R3  E%  AAnorth  +  R4  Er  A Anorth  —  R§  AAsouth  —  RiEzAAsouth  —  RfExAAsouth  =  0  (16) 

One  challenge  to  solving  this  equation  is  how  to  discretize  the  “cross”  terms,  such  as  E ^  and  E%.  The 
approach  used  in  2DFV  consists  of  moving  these  cross  terms  to  the  RHS  of  the  equation  (along  with  the 
electric  field-free  terms)  and  solving  for  them  iteratively.  This  leaves  only  four  electric  field  terms  on  the 
LHS,  corresponding  to  the  normal  field  on  each  of  the  four  faces.  The  resulting  equation  is  as  follows: 


E%  AAeast  ~  Z/W E w A Awest  +  E^ AAnorth  —  R^EftAAsouth  — 

-(Z*  +  Z?E%)AAeast  +  (ZY  +  zYEY)AAwest  -  {R%  +  R?E%)AAnorth  +  (i?f  +  RlEsz)AAsouth  (17) 

The  face-centered  electric  fields  on  the  LHS  of  Eqn.  17  can  be  expanded  using  the  cell-centered  potentials 
as  follows: 
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Cj bi+ 1’j  -  <j>^ 


EEZ  = 
Ef  = 
E%  = 
ESr  = 


Az 

O'--'  -  o'  -'1-'' 

Az 

^i,i+l  _  tpij 

A  r 

A  r 


(18) 


With  appropriate  boundary  conditions,  a  pentadiagonal  linear  system  of  equations  of  the  form  Ax=b 
can  be  constructed  (where  x  is  the  solution  vector  for  the  potential).  This  equation  is  iterated  until  the 
L2-norm  of  the  difference  in  the  potential  field  between  two  iterations  is  less  than  0.001.  For  more  complete 
convergence  studies,  see  Appendix  B.  It  is  important  to  note  that  because  this  potential  solver  is  cell-based, 
the  resulting  electric  field  can  be  calculated  only  on  the  interior  nodes.  Special  treatment  is  necessary  to 
obtain  the  electric  field  on  the  boundary  nodes.  A  flow  diagram  summarizing  2DFV  execution  is  presented 
in  Fig.  2. 


Figure  2.  Flow  Diagram  for  Potential  Solver 


B.  Boundary  conditions  for  linear  system 

Appropriate  boundary  conditions  are  necessary  to  obtain  the  correct  A  matrix  for  the  linear  equation  solve 
step.  Dirichlet  potential  conditions  at  the  anode  and  virtual  cathode  can  be  implement ated  directly  into 
Eqn.  18.  On  the  inner  and  outer  walls,  Eqn.  15  is  modified  by  setting  and  equal  to  zero,  respectively. 
At  corner  points,  both  of  these  conditions  are  satisfied  simultaneously. 

C.  Boundary  conditions  for  electric  field 

As  shown  in  Fig.  2,  separate  boundary  conditions  are  needed  to  establish  the  electric  field  on  the  boundary 
nodes  for  both  normal  and  cross  boundary  conditions.  These  normal  boundary  conditions  correspond  to  Ez 
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conditions  on  the  EW  grid  and  Er  conditions  on  the  NS  grid.  The  Efw  BC  couples  directly  to  the  Dirichlet 
potential  conditions  at  the  anode  and  virtual  cathode.  The  E^s  BC,  representing  the  radial  force  on  the 
plasma  towards  or  away  from  the  walls,  would  be  an  ideal  location  to  implement  a  presheath  correction; 
however,  for  the  results  presented  in  this  paper,  the  E^s  BC  is  simply  the  value  of  Er  on  the  adjacent 
interior  node  of  the  NS  grid. 

The  cross  boundary  conditions  correspond  to  Er  conditions  on  the  EW  grid  and  Ez  conditions  on  the  NS 
grid.  The  E^w  BC  is  calculated  directly  from  the  Dirichlet  potential  conditions  at  the  anode  and  virtual 
cathode.  (For  a  perfectly  conducting  anode,  this  would  translate  to  an  E^w  BC  of  exactly  zero.)  The  E^s 
BC  represents  the  axial  force  on  the  plasma  along  the  dielectric  wall.  Since  no  axial  presheath  physics  are 
considered  in  2DFV,  the  E^s  BC  is  simply  the  value  of  Ez  on  the  adjacent  interior  node  of  the  NS  grid. 

D.  Implementation  details 

The  results  presented  in  this  paper  were  produced  using  MATLAB  7.0.  Sparse  matrix  structures  were  used 
to  store  to  coefficient  matricies  A  and  B.  Solution  time  was  typically  on  the  order  of  5-25  seconds  (for  200- 
750  iterations)  on  a  40x40  node  grid.  When  coupled  into  a  full  hybrid-PIC  code,  it  is  estimated  that  fewer 
iterations  will  be  needed  since  the  initial  guess  for  the  potential  field  from  the  previous  timestep  should  be 
close  to  convergence  if  the  plasma  conditions  are  not  vastly  different.  It  should  also  be  noted  that  MATLAB 
is  a  comparatively  slow  interpreted  language  and  significant  speedup  is  expected  when  this  code  is  ported 
to  Cn — K 


V.  Results  and  Discussion 


Two  major  themes  are  explored  in  this  section  -  the  sensitivity  of  the  solution  to  the  anomalous  mobility 
correction  and  uncertainty  in  Hall  current  calculation.  The  results  presented  in  this  paper  are  based  on  the 
plasma  density,  neutral  density,  ion  velocity,  and  electron  temperature  profiles  presented  in  Fig.  3.  These  are 
instantaneous  plasma  properties  from  a  simulation  of  HPHall,  developed  by  Fife,1  using  a  Bohm- mobility 
correction  to  the  anomalous  mobility.  The  magnetic  field  profile  and  corresponding  electrostatic  potential 
profile  (based  on  the  quasi-ID  potential  solver  in  HPHall)  are  presented  in  Fig.  4. 
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Figure  3.  Initial  data  profiles 
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(a)  Magnetic  field  lines 


(b)  Plasma  potential  (quasi-ID  potential  solver) 


Figure  4.  Field  profiles  from  HPHall 


Classical  mobility 

Traditionally,  the  principal  electron  collision  mechanism  which  impacts  the  electron  mobility  is  electron- 
neutral  elastic  scattering.  The  simple  formulation  used  in  2DFV  is  as  follows: 


v  =  ven  =  2.15 E  —  13  •  nxe 


(19) 


Anomalous  mobility 

The  addition  of  anomalous  mobility  to  Hall  thruster  codes  has  been  attempted  through  both  direct  addition 
to  the  mobility  term  by  Fife:1 


M  M classical  T  Hanom 

and  also  through  increasing  the  effective  electron  elastic  collision  frequency  by  Boeuf:5 


(20) 


P  —  ^en  T  Vanom  (21) 

In  the  framework  of  2DFV,  the  second  type  of  correction  is  used.  The  particular  form  of  this  anomalous 
mobility  correction  is  based  on  the  idea  that  the  electron  can  “jump”  field  lines  by  scattering  off  the  dielectric 
wall  sheath.  The  form  of  the  wall-based  anomalous  mobility  correction  used  in  this  code  is: 

Vanom  =  a-f-  (22) 

ac 

This  form  includes  a  modeling  coefficient  a  and  the  definition  of  the  local  electron  thermal  velocity  as, 


Vth 


8kbTe 

7rrne 


(23) 


A.  Sensitivity 

Potential  profiles  and  mobility  coefficients  for  a  range  of  a  from  0.0  to  0.8  are  presented  in  Fig.  5  and 
Fig.  6.  The  ratio  of  the  parallel  to  perpendicular  mobility  in  Eqn.  9  (assuming  a  ratio  of  fiB  »  1)  can  be 
approximated  by: 


M|| 


(24) 


One  implication  of  this  relationship  is  that  as  the  mobility  ratio  increases  the  potential  solution  should 
closely  resemble  the  magnetic  field  lines.  Classical  mobility  represents  the  lowest  reasonable  bounds  for  v 
and  thus  the  highest  possible  mobility  ratio  expected  in  a  Hall  thruster  simulation.  This  behavior  is  evident 
in  the  resemblance  of  Fig.  5(a)  to  the  magnetic  field  lines  presented  in  Fig.  4. 

An  additional  result  of  note  is  the  change  in  the  orientation  of  the  equipotential  line  (most  evident  in 
Fig.  5(i))  near  the  outer  dielectric  wall  between  0.01  and  0.015  m  from  the  anode.  The  behavior  reflects  the 
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local  drop  in  the  plasma  density  which  is  evident  in  Fig.  3(a).  This  steep  gradient  in  plasma  density  translates 
to  a  strong  Vrp  force  and  a  corresponding  steep  potential  gradient  in  the  radial  direction.  Capturing  such 
2D  effects  is  not  always  possible  with  a  quasi- ID  potential  solver. 

Table  1.  Current  measured  at  Virtual  Cathode 


a 

h  (A 

h  (A) 

h  (A) 

0.800 

1.3057 

4.5503 

5.8559 

0.600 

1.3057 

3.4117 

4.7173 

0.400 

1.3057 

2.2564 

3.562 

0.200 

1.3057 

1.1147 

2.4204 

0.100 

1.3057 

0.63863 

1.9443 

0.075 

1.3057 

0.60419 

1.9098 

0.050 

1.3057 

0.75488 

2.0605 

0.025 

1.3057 

1.7935 

3.0992 

0.000 

1.3057 

19.174 

20.48 

The  anomalous  mobility  correction  has  a  significant  impact  on  the  magnitude  of  the  electron  current  (and 
thus  on  the  discharge  current  and  thruster  efficiency) .  The  sensitivity  of  the  electron  and  discharge  currents 
to  the  anomalous  mobility  correction  are  presented  in  the  Table  1.  The  minimum  in  the  discharge  current 
observed  around  <a=0.075  can  be  interpreted  by  considering  the  following  notional  relationship  between  the 
electron  mobility  and  electric  field: 


Ie  ~  M  •  neE  (25) 

In  the  limit  of  classical  mobility,  there  is  a  significant  potential  drop  in  the  vicinity  of  the  peak  plasma 
density  so  the  neE  term  is  fairly  large.  Furthermore,  the  large  cross  mobility  (n  12)  term  also  contributes  to  a 
very  high  axial  electron  current.  These  effects  diminish  with  the  addition  of  anomalous  mobility.  By  contrast, 
the  addition  of  the  anomalous  mobility  correction  directly  elevates  the  electron  current.  Not  surprisingly, 
somewhere  between  these  limits  for  the  mobility  correction  (in  this  case,  around  <a=0.075),  a  minimum  exists 
in  the  discharge  current.  This  behavior,  while  completely  reasonable  in  the  context  of  a  steady  state  potential 
solver  with  fixed  plasma  conditions,  is  not  necessarily  physically  realistic  in  a  time-dependent  solution  due 
to  the  ability  of  the  plasma  to  respond  to  the  potential  gradient. 


B.  Hall  current 


The  Hall  current  density  is  a  very  difficult  parameter  to  measure  experimentally;  however,  from  simulation 
data  it  can  be  calculated  in  at  least  two  ways.  The  Hall  current  density  can  be  calculated  directly  from 
Eqn.  4  (“jo”)  or  from  the  ExB  drift  as  follows  ( “j Exb” ) • 


j  ExB  =  -qene(E  x  B)\e  =  -qene 


ErBz  -  EzBr 

I  BP 


(26) 


The  integrated  Hall  current  from  the  two  approaches  and  the  corresponding  Hall  parameters  (ratio  of  Hall 
current  to  discharge  current)  are  provided  in  Table  2.  These  results  demonstrate  fundamentally  different 
plasma  behavior.  The  jo  approach  relies  very  heavily  on  the  magnitude  of  the  electron  mobility.  Since  the 
anomalous  correction  directly  raises  or  lowers  this  value,  it  also  serves  to  very  strongly  adjust  the  magnitude 
of  the  j<9-based  Hall  parameter.  The  jExB  approach  does  not  have  this  same  strong  dependence  on  the 
anomalous  mobility  correction.  Since  the  potential  profiles  are  bounded  by  the  anode  and  virtual  cathode 
potentials  and  thus  show  limited  change  with  the  a  parameter,  it  is  not  surprising  that  the  jExB- based  Hall 
parameters  shows  much  smaller  variation  over  the  same  range  of  anomalous  mobility  corrections. 

The  location  of  the  Hall  current  for  the  cases  of  a  =  0.075  and  a  =  0.600,  provided  in  Fig.  7,  clearly 
indicates  very  different  behavior  for  the  two  different  methods  of  calculating  the  Hall  current  density.  When 
calculated  using  jExB ,  the  highest  Hall  current  density  exists  in  the  center  of  the  acceleration  channel  due  to 
the  high  plasma  density  (as  shown  in  Fig.  3(a))  in  the  interior  of  the  thruster.  By  contrast,  when  calculated 
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Figure  5.  2DFV  Results  (Part  1  of  2)  (Note:  Contour  range  not  uniform) 
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Figure  6.  2DFV  Results  (Part  2  of  2)  (Note:  Contour  range  not  uniform) 
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Table  2.  Hall  current  calculations 


a 

Id  (A) 

\Iexb\  (A) 

\Ie\  (A) 

^ ExB 

^ 0 

0.800 

5.86 

8.33 

15.38 

1.4 

2.6 

0.600 

4.72 

8.87 

16.68 

1.9 

3.5 

0.400 

3.56 

9.73 

19.17 

2.7 

5.4 

0.200 

2.42 

10.92 

26.69 

4.5 

11.0 

0.100 

1.94 

11.43 

43.08 

5.9 

22.2 

0.075 

1.91 

11.50 

54.92 

6.0 

28.8 

0.050 

2.06 

11.52 

80.83 

5.6 

39.2 

0.025 

3.10 

11.52 

179.48 

3.7 

57.9 

0.000 

20.48 

13.19 

14404.00 

0.6 

703.3 

using  jo,  the  Hall  current  density  (both  positive  and  negative)  is  preferentially  concentrated  near  the  inner 
and  outer  dielectric  walls.  The  large  magnitude  of  the  Hall  current  density  in  these  regions  is  largely  due  to 
the  high  radial  current  densities  near  the  walls. 


(a)  j ExB  ,  ol  =  0.075 


(b)  jd,a= 0.075 


(c)  j  ExB,  O'  =  0.600 


(d)  jo,  ol  —  0.600 


Figure  7.  Hall  current  density  profiles  (Note:  Contour  range  not  uniform) 


The  significant  differences  in  Hall  current  density  behavior  displayed  by  these  results  must  be  put  into 
context.  To  begin,  this  particular  plasma  configuration  (and  the  resulting  Hall  current  profiles)  are  based 
on  computational,  rather  than  experimental,  initial  conditions.  Furthermore,  although  this  is  a  steady-state 
potential  solver,  the  ability  of  the  plasma  to  react  to  the  electric  field  implies  that  this  potential  state 
is  not  necessarily  a  true  steady-state  configuration  that  might  be  observed  experimentally.  Even  taking 
these  considerations  into  account,  some  mechanism  to  account  for  anomalous  transport  without  drastically 
modifying  the  electron  mobility  coefficient  might  help  to  reduce  the  discrepancy  between  the  Hall  current 
density  results  presented  in  this  paper.  An  possible  extension  of  this  2D  axisymmetric  formulation  to  account 
for  j  x  B  plasma  turbulence  effects  directly  is  derived  in  Appendix  A. 

Further  issues  with  the  calculation  of  the  Hall  current  from  the  j  ExB  drift  alone  appear  when  considering 
an  additional  component  to  the  azimuthal  electron  drift  from  the  diamagnetic  drift  of  the  electrons: 
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(27) 


Vp  x  B 

VDe  qeneB2 

The  calculated  ExB  and  diamagnetic  drifts  for  an  a  =  0.6  simulation  are  provided  in  Fig.  8.  The 
integrated  ExB  electron  current,  as  shown  in  Table  2,  is  8.87  A  while  the  integrated  diamagnetic  electron 
current  is  -4.66  A.  Close  inspection  of  Fig.  8  shows  that  the  diamagnetic  current  density  peaks  (negative)  near 
the  anode  and  while  the  ExB  current  density  peaks  (positive)  much  closer  to  the  exit  of  the  thruster.  Again, 
in  interpreting  these  results,  it  is  crucial  to  understand  that  all  results  presented  in  this  paper  are  based 
on  a  single  set  of  computationally  derived  input  data.  For  this  case  in  particular,  the  very  strong  negative 
diamagnetic  electron  drift  is  likely  due  to  the  strong  gradient  in  electron  temperature  near  the  anode  seen 
Fig.  3(d).  While  the  inclusion  of  the  diamagnetic  electron  drift  does  little  to  explain  the  discrepancy  between 
the  jo  and  JexB  drifts,  the  magnitude  of  the  diamagnetic  electron  drift,  even  for  this  “non-experimental” 
case,  indicates  that  it  could  be  a  significant  drift  mechanism  in  realistic  discharge  physics. 
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Figure  8.  Classical  fluid  drifts  perpendicular  to  B 


VI.  Conclusion 

This  work  demonstrates  the  formulation  of  a  finite-volume  discretization  scheme  to  solve  a  2D  axisym- 
metric  version  of  the  current  conservation  equation.  The  resulting  simulation  code,  2DFV,  is  used  to  model 
the  acceleration  channel  of  a  Hall  thruster.  Performance  of  this  code  indicates  that  it  is  well  suited  to 
update  existing  quasi- ID  potential  solvers  in  hybrid  PIC  simulations  (although  a  corresponding  2D  electron 
energy  solver  is  also  desirable  for  true  2D  capability).  2DFV  results  presented  in  this  work  display  physically 
reasonable  behavior  in  the  limit  of  classical  mobility  and  demonstrate  a  capacity  to  incorporate  anomalous 
mobility  corrections  into  the  computational  framework.  Results  provided  indicate  that  strong  discrepancies 
exist  between  the  azimuthal  electron  current  as  calculated  from  the  2DFV  framework  and  from  the  azimuthal 
component  of  the  ExB  drift.  Finally,  an  extension  of  the  2D  axisymmetric  formulation  to  account  for  j  x  B 
plasma  turbulence  effects  is  proposed  in  Appendix  A. 

Appendix  A  -  Formulation  with  plasma  turbulence  terms 

Using  basic  concepts  from  linear  perturbation  theory,  it  is  possible  to  include  “plasma  turbulence”  terms 
into  a  Generalized  Ohm’s  law  formulation  by  replacing  individual  field  variables  with  a  steady-state  term 
plus  a  fluctuation  term  as  follows  (removing  vector  notation  where  applicable  for  simplicity): 


je 

=  Je  +  j'e 

B 

=  B  +  B' 

E 

=  E  +  E' 

V 

=  p  +  p' 

(28) 

Plug  these  into  the  Eqn.  2  and  taking  the  time  average  (note  that  the  time  average  of  the  fluctuating 
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components  is  zero,  so  the  electric  field  and  pressure  terms  contain  only  steady  state  components)  leaves  the 
following  formulation  for  the  electron  current  density: 

je  =  Mie  X  B)  +  M  <  j'e  x  B'  >  EpneE  +  pVp  (29) 

The  additional  cross  correlation  term,  <  j'e  x  B'  >,  represents  the  time  average  of  the  cross  product 
interaction  between  oscillations  in  the  electron  current  density  and  the  magnetic  field.  Depending  on  the 
particular  model  used  to  represent  this  term,  it  can  provide  a  spatially  varying  positive,  zero,  or  negative 
contribution  to  the  electron  current  density.  Modeling  this  term  is  beyond  the  scope  of  this  effort;  however, 
future  research  models  for  Bohm  mobility  should  couple  much  more  naturally  with  this  approach  than  with 
2DFV  due  to  the  explicit  inclusion  of  these  turbulent  transport  terms.  To  simplify  the  equation  system,  the 
<  j'e  x  B'  >  term  will  be  represented  by  the  vector,  (az,  <ar,  a#),  corresponding  to  the  turbulent  transport 
corrections  in  the  three  principal  axes.  Expanding  the  full  electron  current  density  equation,  a  set  of  three 
equations  can  be  written  for  the  electron  current  in  the  axial,  radial,  and  azimuthal  directions  (dropping  the 
subscript  e  for  convenience). 

jz  =  p{jeBr  jrBo  +  olz )  -b  pucEz  T  /iV zp 

jr  nUzBo  jeBz  T  C^r)  T  fATlcEr  T  fl W rp 

je  =  vtirBz  -  jzBr  +  OLe)  +  pneEe  +  pVop  (30) 

Due  to  the  axisymmetric  geometry  of  the  Hall  thruster,  the  electric  field  and  gradient  of  the  pressure 

in  azimuthal  direction  are  identically  zero.  With  this  simplification,  the  azimuthal  current  is  no  longer  a 
function  of  any  azimuthal  quantities  and  is  written  as: 


je  =  n(jrBz  -  jzBr  +  ae) 


(31) 


Now  it  is  possible  to  substitute  jo  into  the  equations  for  the  axial  and  radial  electron  current  density. 
This  leads  to  the  following  set  of  coupled  linear  equations: 


jz  =  p{p{jrBz  -  jzBr  +  ae)Br  -  jrBe  +  az)  +  pneEz  +  /xV zp 

jr  =  n(jzBo  p^jrBz  jzBr  ae)Bz  T  ov)  T  iiTicEr  T  /rVrp  (3^) 

This  equation  system  can  be  solved  to  isolate  the  axial  and  radial  components  of  the  electric  current  density: 


jz  =  Mil  {li2aeBr  +  [iolz  +  pneEz  +  pSJ  zp) 

+pi2(fJ2aoBz  +  (iOLr  +  fmeEr  +  pSJ  rp) 
jr  =  M21  {pi2  olq  Br  +  piOLz  +  pnfieEz  +  p'VzP) 

TM22  {pi2OL$Bz  +  pOLr  +  pineEr  +  pEJ rp)  (33) 

Where  the  mobility  coefficients  are  defined  in  Eqn.  8  and  Eqn.  7  and  the  form  for  the  azimuthal  electron 
current  is  defined  in  Eqn.  31. 


Appendix  B  -  Convergence 

Different  L 2  norm  convergence  criteria  are  applied  to  a  40x40  node  simulation  with  a  =  0.6.  The 
resulting  currents  are  provided  in  Table  3.  To  assess  grid  convergence,  internal  plasma  quantities  (electron 
temperature,  ion  velocity,  plasma  density,  and  neutral  density)  are  linearly  interpolated  to  the  desired  grid 
for  a  simulation  with  an  L2  norm  convergence  criteria  of  0.001  and  an  anomalous  correction  of  a  =  0.6.  This 
use  of  interpolated  plasma  properties  can  be  anticipated  to  perturb  grid  convergence;  however,  as  the  results 
in  Table  4  demonstrate,  approximate  grid  convergence  can  be  achieved  with  a  fine  enough  mesh.  (Note:  the 
20x20  grid  is  very  accurate  because  it  does  not  require  linear  interpolation  from  the  40x40  grid.) 
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Table  3.  L2  norm  convergence  criteria 


L2  norm 

Ie  (A) 

Id  (A) 

Iterations 

Time  (sec) 

10 

3.4336 

4.7393 

200 

6.1 

1 

3.4138 

4.7195 

316 

9.7 

0.1 

3.4119 

4.7175 

430 

13.1 

0.01 

3.4117 

4.7173 

546 

16.3 

0.001 

3.4117 

4.7173 

660 

19.8 

0.0001 

3.4117 

4.7173 

774 

23.4 

Table  4.  Grid  convergence  criteria 


Grid  size 

Ie  (A) 

Id  (A) 

10x10 

4.0189 

5.3135 

20x20 

3.411 

4.7156 

25x25 

3.4076 

4.712 

30x30 

3.3825 

4.6872 

35x35 

3.3859 

4.6909 

40x40 

3.4117 

4.7173 
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